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Summary 


In this paper, we propose an optimal least-squares finite element method for 2D and 
3D elliptic problems and discuss its advantages over the mixed Galerkin method and 
the usual least-squares finite element method. In the usual least-squares finite element 
method, the second-order equation —V • (Vu) + u = / is recast as a first-order system 
— V • p + u = /, Vu — p = 0. Our error analysis and numerical experiments show 
that, in this usual least-squares finite element method, the rate of convergence for flux 
p is one-order lower than optimal. In order to get an optimal least-squares method, the 
irrotationality V x p = 0 should be included in the first-order system. 
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1 . Introduction 


The least-squares finite element method (LSFEM) discussed here is based on minimizing 
the norm of the residuals of partial differential equations. In order to use elements, 
the second-order 2D elliptic partial differential equation is reduced to a system of three 
first-order differential equations by introducing two more unknown variables (flux). This 
idea was first proposed by Lynn and Arya[l] and Zienkiewicz[ 2 ], and was an important 
contribution to the development of least-squares finite element methods. This procedure 
has long been considered as a standard way to develop least-squares methods. 

In this paper, we present both theoretical analysis and numerical results to show that, 
this simple procedure of reduction destroys ellipticity and the usual LSFEM is not optimal, 
that is, the rate of convergence for flux is one-order lower than optimal. In order to get 
an optimal LSFEM, the compatibility condition (the irrotationality) should be included in 
the first-order system. 

The plan of the presentation is as follows. In Section 2 , we introduce the model 
problem and notations. In Section 3, we give a short summary on the related mixed 
Galerkin method for the purpose of comparison. In Section 4 , we analyse the usual LSFEM 
and explain where the trouble comes from. In Section 5 , an optimal LSFEM and the error 
estimation are presented. In Section 6 , we discuss how to deal with some more general 
elliptic problems. In Section 7, numerical results are given. 


2. Preliminaries and Notations 

In this paper, we present the essential idea of the optimal least-squares finite element 
method by solving the following second-order elliptic boundary-value problem: 

—V • Vu + u = /(x) in fi, 

Vix-n=^(x) on T, ^ 

where fl C S" (n = 2 or 3) is an open bounded convex domain with a piecewise C 1 
boundary T, x = (%i ,x 2 ,X3) is a point in fi, n = ( 111 , 712 , 713 ) is a unit outward normal 
vector on the boundary, and /(x) and g(x) are given functions. Without loss of generality, 
we shall hereafter consider only the homogeneous boundary condition for simplicity, that 
is, we shall take <?(x) = 0. The primal variable u can be, for instance, temperature for 
heat conduction; potential for incompressible and irrotational flow; or electric potential for 
electromagnetics, etc. 

Throughout this paper, we use the following notations. L 2 (Q) denotes the space of 
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square-integrable functions defined on Q equipped with the inner product 



uvdx 


u,v £ Z 2 (^) 


and the norm 



= («»“) 


u £ £12(0). 


denotes the Sobolev space of functions with square-integrable derivatives of order 
up to r. | • | r and || • || r denote the usual seminorm and norm for H r (Cl) y respectively. For 
vector-valued function p with n components, we have the product spaces 


(£»(«))"> (fl'fs))’, 


and the corresponding norm 


liti-Elnli;. IIpIIJ = E llPill?- 

3 = 1 3 = 1 

Further we define the function spaces 

H = {v £ 

5={qe (^ 1 (n))"|q ■ n = 0 on T}, 

W = {q € (i 2 (ft))"|q • n = 0 an T}, 

and the corresponding finite element subspaces Hk, Sk and Wh, i-e., Hh, and Sh are the 
spaces of continuous piecewise polynomial functions of order k, and Wh is the space of 
continuous piecewise polynomial functions of order k — 1. Here the parameter h represents 
the maximal diameter of the elements. By the finite element interpolation theory [3, 4], 
we have: Given a function u £ and a function p £ there exist an 

interpolant u h £ Hh and p £ Sh such that 

||ii-«‘||o<Cfc‘ +, HU +1 , 

||«-i‘lll <Ck‘|MU + i, (2) 

||p-p‘|lo<Cfc‘ +1 ||p|| k+1 , 

|[p-p h |l> <ct**IIp|Im-,, 

here and below C denotes a constant independent of the mesh parameter h y with possibly 
different values in each appearance. 
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We would also like to write down Green’s formula 


(V -q,v) + (q, Vv) = L uq • nds. 


( 3 ) 


3. The Mixed Galerkin Method 

The most commonly used method for problem (1) is the classical Galerkin method. 
However, a posteriori numerical differentiation is required to obtain dual variables (flux 
for heat transfer; velocity for fluid flow; or electric field intensity for electromagnetics) 
which are often of most interest. In general, the accuracy of so computed dual variables is 
one-order lower than that of the primal variable. 

Mixed Galerkin methods were devised in the hope of getting better accuracy for dual 
variables [3, 5]. Here the term “mixed” means that both the primal variable and the dual 
variables are approximated as fundamental unknowns. 

In mixed methods, problem (1) is decomposed into an equivalent first-order system: 


-V • p + u-f 

in f), 



Vu — p = 0 

in ft, 


(4) 

p • n = 0 

on r. 



Then the Galerkin principle is applied to problem 
weak statement: Find (u,p) G H xW such that 

(4). This 

leads to the mixed Galerkin 

1 ( uv + Vv • p)dx = / fvdx. 
J n Jn 

/ (Vu • q — p • q)dx = 0 
Jn 

in ft 
in ft 

Vug h, 
Vq G W. 

(5) 

It is well known that problem (5) corresponds to a saddle-point variational problem, 
and thus in order to guarantee the existence of the solution, the pair (u,p) must satisfy 
the following Babuska-Brezzi condition: 

ueH ( / ^-P^XIMli) 1 > C'HpHo 

u#0 , ' n 

Vp G W. 

(6) 
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The Babuska-Brezzi condition precludes the application of simple equal-order finite 
elements. It can be proved that the finite element spaces H h and Wh satisfy the dis- 
crete Babuska-Brezzi condition (6), and if the solution (u,p) of (4) belongs to H k+1 (Sl) x 
(ff k (fl)) n , we have the following error estimate[5|: 

l« - “/tli + Up - Pfcllo < ch k (j u j k+1 + ||p|| fc ). (7) 

The estimate (7) tells us that in this mixed method the accuracy for flux p is always 
one-order lower than that for the primal variable u. 

By inspecting equation (5), we know that the matrix associated with the mixed method 
is non-positive. This makes the use of iterative methods to solve large-scale problems very 
difficult. 


4. The Usual LSFEM 

The usual LSFEM [1,2] is also based on first-order system (4). For 2D problems, the 
first-order system in (4) consists of three equations and three unknown functions. The first- 
order system with an odd number of unknowns and an odd number of equations cannot 
form an elliptic system in the ordinary sense. For 3D problems, although the first-order 
system in (4) has four unknowns and four equations, it is easy to verify that the system is 
not elliptic in the ordinary sense. This fact makes us suspect that the LSFEM based on 
system (4) will not be optimal. 

Now let us analyse the usual LSFEM which minimizes the following functional: 

I : H x S — R, 

T ( u > P) = II - V ■ p + u - /||2 -(- ||Vu - p|||. (8) 

Taking variation of I with respect to u and p, and letting 81 = 0, 8u = v and ^p = q lead 
to a least-squares weak statement: Find U - (u,p) G H x S, such that 

B(U,V)=L(V) V7 = ( M )eFx5, (9) 


where 

V) = (~ V • p + «, -V ■ q + v) + (Vu - p, Vv — q), 

uy) = (/, -V • q + v). 
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The corresponding finite element problem is then to find Uh — (u^, Ph) E Hh x S&, 
such that 

*(ff k ,V fc ) = L(V h ) VV h = (t*,q fc ) eH h x S h , (10) 


where 

B(Uh,Vh ) = (-V • pa + U ) i,— V • + ua) + ( Vti a - Pa, Vua - q^), 

L(VA) = (/,-V-qA+ V A). 


It is easy to verify that 

B(C,V)<C||P|M|V||, (11) 

where ||{7|| 2 — ||u|| 2 + IIpIIo + ||V ■ p||g. Thus, B(U, V) is continuous on H x S. Since 
B(U,V) is symmetric, the inequality in the Lax-Milgram theorem reduces to the single 
coercivity requirement: There exists a constant a > 0 such that for V G H x S 

B(V,V)>a ||Vf. (12) 

Let us prove the coercivity (12). We know that 

B(V, V) = II - V ■ q + HI? + ||V» - q||o- (13) 

Consequently, 

B(V,V) > ||V ■ q - I>||; = ||V ■ q(|o + ||v|i: - 2(V ■ q,„), (14) 

B(V, V) > ||V» - q|| o = l|Vt»||S + ||q||o - 2(q, Vr). (15) 

By combining (14) and (15) together and using Green’s formula (3) and the boundary 
condition, we obtain the coercivity: 

B(V, V) > 0.5(|| Vv||S + IHIJ + IIV- q||J + ||q||J), 

or 

B(v,v)>oj(iHiJ + ||v.q|i; + ||q|i;). (is) 

REMARK. We may say that the coercivity (16) is incomplete or deficient, because the 
derivatives of q are not completely controlled in general. This explains why the usual 
LSFEM often does not behave well. 

Therefore, the following theorem about the rate of convergence of the usual least- 
squares finite element solutions can be derived. 
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THEOREM 4-1- Assume that /(x) £ L 2 (H), the solution (u,p) of (4) belongs to ,ff fc+1 (D)x 
(-ff fe+1 (0)) n , and the finite element interpolation estimates (2) hold. Then for the approx- 
imate solution associated with (10), we have the error estimate: 

11“ - “/till + ||V • (p - Pfc)||o + Up - Philo < C^flMU+i + ||p||fc+i). (17) 

This theorem shows that the accuracy of the flux p for the usual LSFEM is one- 
order lower than optimal. Even so, the usual LSFEM has significant advantages over the 
mixed Galerkin method. Namely the usual LSFEM is not subject to the Babuska-Brezzi 
condition and thus can accommodate simple equal-order elements, and the resulting matrix 
is symmetric and positive definite and thus simple iterative methods, such as conjugate 
gradient methods, can be employed and vectorization and parallelization are trivial. 


5. The Optimal LSFEM 

Conservative laws and constitutive laws in physics are in general governed by first- 
order systems. For historic reasons (convenience for hand calculation and analysis), the 
equations in a first-order system are combined into a high-order partial differential equation 
(or equations) with one or less unknowns. For example, for incompressible and irrotational 
flows, by introducing the potential (or the stream function in 2D cases), the incompressibil- 
ity and the irrotationality are combined into a second-order Laplace or Poisson equation. 

We believe that in the computer age this transformation is unnecessary. We may 
solve directly the original first-order governing equations by LSFEM. For flows considered 


in this paper, the governing equations are the following: 

— V-p + u = / in 0, (18-1) 

V x p = 0 in Cl, (18.2) 

Vu — p — 0 in Cl, (18.3) 

p ■ n = 0 on T. (18.4) 


Here (18.1) is the mass conservation, (18.2) represents the irrotationality, and (18.3) is the 
constitutive relation. 

For 2D problems, the first-order system in (18) consists of three unknowns and four 
equations. At first glance, one may think that this is an overdeter mined system. In fact, 
this is a determined system in the sense that the components p\ and P 2 of p in (18.3) are 
not completely independent. They must satisfy (18.2). 
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Let us show that system (18) is determined and elliptic. In 2D cases, we introduce 
a dummy variable 0 (this technique was first pointed out to us by C.L. Chang for 2D 
problems), and rewrite system (18) as 


dx\ 

+ 

dx 2 

— 

dpi < 

dx 2 ( 

du 

d(f> 

dx\ 

dx 2 

du 

Ck 

d<f> 

* Ck 


-u = -f 
^ = 0 


- Pi - o 


-p 2 = 0 


+ P2^2 — 0 

0 = 0 on r. 


in ft, 

in ft, 
in II, 
in ft, 

on r. 


(19.1) 

(19.2) 

(19.3) 

(19.4) 

(19.5) 

(19.6) 


5(19.3)/$x 2 — $(19.4)/5xi leads to d 2 <f>/dx 2 + d 2 <f>/dx\ = 0. We have already specified 
that 0 = 0 on T, thus 0 = 0 in 0. This means that system (18) with three unknowns and 
four equations is indeed equivalent to system (19) with four unknowns and four equations. 

Now we write (19) in a standard matrix form: 


in which 


Since 


Ai = 


A = 


. #U dll . 

Ai bA 2 - (-Au^f in ft, 

OX i OX2 


0 10 0 

-10 0 0 
0 0 0 -1 
0 0 10 


u 


/I 

0 

0 

°\ , 

1 0 

1 

0 

0 I . 

° 

0 

1 

0 i ’ A 2 - I 

\0 

0 

0 

1/ ' 


/ 0 0-10 

0 0 0 0 

-10 0 0 
\ 0 -1 0 0, 


f = 



det(Aj£ + A 2 t/) = dei 


0 

0 



= (c + v r t o 


( 20 ) 


for all nonzero real pair (£,?/), system (19) and thus system (18) is determined and elliptic, 
as contended. 


8 



In 3D cases, by introducing dummy unknowns 0, uii , u? 2 , , we write the following 

first-order system with eight unknowns and eight equations: 


in SI, 
in fl, 
in 0, 
in fi, 
in II, 
in Si, 
in Q, 


dp 1 

+ ■ 

dp 2 

. dpz 



-/ 

dx\ 

dx 2 

+ • 
( 

dx 3 

— 11 = 

dp 2 

+ 

dp 3 

+ 

00 

+ 


= 0 

dx 3 

dx 2 

0a; 1 

= 

dp 3 

+ 

dpi 

+ 

00 

+ 


= 0 

dx 1 

dx 3 

dx 2 

W 2 = 

dp 1 

+ 

dp 2 

+ 

00 

+ 


= 0 

dx 2 

dx\ 

dx 3 

0^3 = 

du> 2 


d U> 3 


du 




dx 3 

+ 

dx 2 

+ 

dx 1 

— 

pi : 

= 0 

du> 3 


du>i 


du 




dx 1 

+ 

dx 3 

+ 

dx 2 

— 

P2 = 

= 0 

dw 1 


du> 2 


du 




dx 2 

+ 

dx 1 

+ 

dx 3 

— 

P3 = 

= 0 


du >2 du> 3 — 
dxi dx 2 dx 3 
PiTii + p 2 n 2 +pzn 3 — 0 
0 = 0 

U>1 = U>2 = = 0 


in O , 

on r, 
on r, 
on r, 


( 21 ) 


It is not difficult to verify that 0 — — u) 2 — u) 3 — 0, and thus system (18) with four 

unknowns and seven equations is indeed equivalent to system (21) with eight unknowns 
and eight equations. 


We may write system (21) in a standard matrix form: 


in which 



+ A 2 


0u 

dx 2 


+ A 3 


0u 

dx 3 


+ Au = f 


in SI, 


( 22 ) 


Ai 


/I 0 0 

0 0 0 
0 0-1 
0 10 
0 0 0 
0 0 0 
0 0 0 
\0 0 0 


0 0 0 0 
0 10 0 
0 0 0 0 
0 0 0 0 
1000 
0 0 0 0 
0 0 0 1 
0010 


0) 



0 / 


/ 0 1000 0 0 0\ 

0 0100 0 00 ' 

0 0001 0 00 
-10000 0 00 
0 0000 0 01 ’ 
0 0010 0 00 
0 0000-100 
Vo 0000 0 10/ 
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/0 

0 

1 

0 

0 

0 

0 

0\ 


0 

0 

0 

-1 

0 

0 

0 

°\ 

0 

-1 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

A — 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

-1 

0 

> A — 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 


0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 j 


0 

0 

-1 

0 

0 

0 

0 

0 

\o 

0 

0 

0 

0 

0 

0 

1 ) 

1 

\ 0 

0 

0 

0 

0 

0 

0 

0 / 


Since 


f = 


/ * 

0 

C 


<fef(A!£+A 277 +A 3 £) = det 


0 
0 
\ 0 




/Pl\ 

0 1 


P 2 

0 


P 3 

0 


u 

0 

» u = 

<!> 

0 


U >1 

0 

\ r\ i 


u >2 

l J 


V 

-C 

0 


-v £ 

0 0 


C 

V 

-£ 0 
0 0 


0 0 

0 £ 


0 

0 

0 


V 

c 

0 £ 0 

0 T] 0 

0 £ 0 


.Wj . 


0 

0 

0 

0 

0 

c 


0 0 0 £ 


0 

0 

0 

0 

-c 

0 

£ 

V 


0 \ 

0 

0 

0 

V 

-£ 

0 

u 


= -(£ +*? +< ) * 0 


for all nonzero real triplet (£, 77 ,^), system ( 21 ) and thus system (18) is determined and 
elliptic. 


Now we may use the error analysis developed by Aziz, Kellogg and Stephens[ 6 ] for 
general elliptic systems to show that the LSFEM based on system (18) is optimal. However, 
their method involves high-level mathematics. Here we use elementary analysis to give a 
poof of the optimality. The key point of this proof is the following technical lemma: 

LEMMA 5.1. Every function q £ S = {q G |q • n = 0 on T} satisfies: 

N?<l|V-q||5 + ||Vxq||J. (23) 


The lemma is discussed in Girault and Raviart[7], and the complete proof can be 
found in Grisvard[ 8 ]. Here we give an elementary proof for 2D rectangular domains. For 
3 D rectangular domains, the proof is similar. 


Proof. 
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Since 


l|V • q||o + IIV X q||g = |q|f + f(pL^ + pLpl-pLpL_ fafSL)*, 

Ju uxi ox 2 ox i a*2 dxi dx 2 dx\ dx 2 ' 

using integration by parts, we have 

l|V-q||S + ||Vxq||2 = |q|;+ /( ni „ + „ 2{2 |?L _ „ l9J |?1 . . 

,/r Ox 2 axj (7x2 


n 2 ?i 


C?Xl 


)ds. 


The boundary term in the foregoing equation is equal to zero by virtue of the boundary 
condition q ■ n = 0. For example, for the part of boundary with n j =1, n 2 = 0, the 
boundary condition is q\ = 0, and thus dq\/dx 2 = 0. Therefore, the boundary integration 
associated with this part of boundary is equal to zero. 


The optimal LSFEM minimizes the following functional: 

/ : H x S -+ ft, 

J (“> P) = || - V • p + u - f\\l + || V x p||j} + ||Vu - p||o - (24) 

Taking variation of I with respect to u and p, and letting SI — 0, Su — v and ^p = q lead 
to a least-squares weak statement: Find U = (u, p) e H x S, such that 

B(U, V) = L(V) VF = {v,q)eHx S , (25) 

where 

B(U, F) = (- V ■ p + u, -V • q + u) -f- (V x p, V x q) + (Vu - p, Vu - q), 

= q + v). 

The corresponding finite element problem is then to find Uh = (uh,Ph) 6 H h x Sh, 
such that 

B(U h ,V h ) = L(V h ) WV h = (v h ,q h ) <= H k x S h , (26) 

where 

B(Uh,Vh ) = (—V • p h + Uh, —V • q^ + Vh ) + (V x Ph, V x q/,) + (V«a — p^, Vvh — q^), 

L(V h ) = (f,-V- q h + v h ). 

It is easy to verify that 

B(U,V)<C\\U\\-\\V\\, (27) 
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where 


IW = INI? + IIpII? 


Thus, B(U, V) is continuous on H x S. Since B(U,V ) is symmetric, the inequality in the 
Lax-Milgram theorem reduces to the single coercivity requirement: There exists a constant 
a > 0 such that for V 6 H x S 

B(V,V)>a IIVII 2 . (28) 

Following the similar argument as in Section 3, we may get 

B(V, V ) > 0.5(11 Vvj|g + ||o||: + ||V • q||; + II V x q||> + ||q||’). (29) 

The combination of (29) and Lemma 5.1 yields the coercivity (28). 

Once the coercivity is proved, the derivation of the following theorem is trivial. 

THEOREM 5.2. Assume that /(x) £ L 2 (ft), the solution (it, p) £ iT fc+ 1 (fi) x (ff fe+1 (fl)) n 
and the finite element interpolation estimates (2) hold. Then for the approximate solution 
associated with (26), we have the error estimate: 

||« - “fclli + ||p — p/i||i < C'h fc (||u||fc +1 + !|p|| fc-pi ). (30) 

This theorem shows that the rate of convergence (in H 1 norm) of the LSFEM based 
on the full first-order system (18) is optimal for all variables. The optimal L 2 convergence 
can be obtained by Aubin-Nitsche method[3j. The optimality attributes to the fact that 
the optimal LSFEM controls not only the divergence, but also the curl of the error of flux. 


6. Discussion 


If there is no u term in the first equation of (1), the corresponding (4) becomes the 
so called div-grad system and the corresponding (18) is the div-grad-curl system. If the 
boundary condition is still only related to flux p, then the calculation of u and p can be 
separated. We may use the LSFEM based on the div-curl system to obtain p first, then 
use p to calculate u. The LSFEM based on the div-curl system is optimal [9, 10]. 

We also would like to discuss the optimal LSFEM for some more general cases. For ex- 
ample, 2D seepage can be modelled by a system of first-order partial differential equations 
of the form: 

dp i 


dxi 


+ 


Pi — ®ii 


dPi f - Q 

= / m SI, 

Ox 2 

(31.1) 

u du . 


+ a i2 q m D, 

'1 OX 2 

(31.2) 
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(31.3) 


P 2 = «21 


du 

dx\ 


+ 022 


du 

dx2 


in Q, 


PiTii +p 2 n 2 = g on T, (31.4) 

where u denotes the hydraulic head, p 2 and p 2 are the components of seepage velocity, and 
f and g are given functions. 

For simplicity, we consider the case in which the coefficient a^j are constant, and the 
matrix (a^) is symmetric and positive definite. Thus, the constitutive relations (31.2) and 
(31.3) are invertible: 


du 
dx 1 


22 — o,j 2 a 21 ) 


( a 22 Pi ~ O-I 2 P 2 ), 


du 1 / 

^2 = (ana J2 -a 12 a 21 ) (flllP2 ~ C2lPl) - 

From (32.1) and (32.2) we can obtain the compatibility condition: 


a 22 


dpi 
dx 2 


— a 12 


dp 2 
dx 2 


dp 2 dpi 

~~ a ll a H a 21 Q — 0. 

dx 1 


(32.1) 

(32.2) 


(33) 


The LSFEM based on equation (31.1), (31.2), (31.3) and (33) will be optimal. 3D 
seepage can be similarly treated. 


7. Numerical Results 


As the first example, we chose 


/-(2.-ix£-£)+(a»-iXT-T> + (£-xM£-£ ) -n, 


where SI — {( x,y ) 6 -ft 2 : 0 < x < 1, 0 < y < 1} is the unit square with the boundary T. 
The boundary conditions are 


Pi - 0 on r 2 = {(x,y) £ r : x = 0}, 

pi = 0 on r 3 = {( x,y ) £ r : x = 1}, 

P2 = 0 on r 2 = {(a:, p) £ T : p = 0}, 

P2-0 onT 4 -{(x,y)£ r-.y = I}. 
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The exact solution should be 


_ (— _ -mO- _ 

U 2 3 2 3 ’ 

^ = (l/ 2 -y)(y-y). 

Numerical experiments were carried out using bilinear elements on uniform meshes with 
1/h — 4,9,20,29. We calculated the L 2 errors for u and p: 

e u = ||U -ttfcHo, e p = (||pi ~Pih\\l + ||P 2 -P2fc||o)^- 

The numerical results on the rates of convergence are given in Fig.l (a) and (b). As 
expected, the rate of convergence of flux p for the usual LSFEM is 0(h) } which is one- 
order lower than optimal. The rates of convergence are 0(h 2 ) for both the primal variable 
u and the dual variables p for the optimal LSFEM. 

As the second example, we chose 

f = (57t 2 + l)cos(27rz)co.s(7ry) 

with the same boundary conditions as in example 1. The exact solution is 

u = —cos(2'Kx)cos{Try\ 

p 1 = 27T3m(27rx)cos(7ry), 
p 2 — 7rcos(27rx)sin(7ry ) . 

The numerical rates of convergence are shown in Fig.l. (c) and (d). It seems that the 
rate of convergence of p for the usual LSFEM is getting better when the mesh is refined. 
However, the error of p for the usual LSFEM is quite large. 

Here we should mention that in all of our calculations, 2x2 Gaussian quadrature 
was used for finite element solutions, and 3x3 Gaussian quadrature was used for error 
evaluation. The LSFEM with numerical quadrature is equivalent to a weighted colloca- 
tion least-squares method. We may use this idea to choose an appropriate number of 
Gaussian points. The usual LSFEM with one-point quadrature will not work, because it 
corresponds to solving a underdetermined algebraic system. The optimal LSFEM with 


14 


one-point quadrature works. However, the computed nodal values of p have oscillations; 
but the values of p at Gaussian points are correct. 
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